The gradient of Antarctic Ice Sheet
Contents
The gradient of Antarctic Ice Sheet¶
As an aside, in this section we use the concept of gradient introduced on the previous page to compute the slope of the Antarctic Ice SHeet using some surface elevation data.
In this section we will:
loas the bedmachine dataset from a version stored in the cloud
computes the surface gradient in the \(x\) and \(y\) directions
computes the slope from these two components
Import packages¶
We will import packages required for opening the data (fsspec), doing computations (numpy), storing the data and results in memory and plotting them (xarray), and plotting with interactive features (hvplot.xarray).
import fsspec
import numpy as np
import xarray as xr
import hvplot.xarray
Lazily load the data.¶
First we get a ‘mapper’ which points to the data in cloud storage, then we lazily load it in to a multi-dimensional array called an xarray, which we cam bm.
bm_mapper = fsspec.get_mapper('gs://ldeo-glaciology/bedmachine/bm.zarr', mode='ab')
bm = xr.open_zarr(bm_mapper)
Let’s take a look at the xarray we get from the cell above
bm
<xarray.Dataset>
Dimensions: (y: 13333, x: 13333)
Coordinates:
* x (x) int32 -3333000 -3332500 -3332000 ... 3332000 3332500 3333000
* y (y) int32 3333000 3332500 3332000 ... -3332000 -3332500 -3333000
Data variables:
bed (y, x) float32 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
errbed (y, x) float32 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
firn (y, x) float32 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
geoid (y, x) int16 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
mask (y, x) int8 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
source (y, x) int8 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
surface (y, x) float32 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
thickness (y, x) float32 dask.array<chunksize=(3000, 3000), meta=np.ndarray>
Attributes: (12/25)
Author: Mathieu Morlighem
Conventions: CF-1.7
Data_citation: Morlighem M. et al., (2019), Deep...
Notes: Data processed at the Department ...
Projection: Polar Stereographic South (71S,0E)
Title: BedMachine Antarctica
... ...
spacing: [500]
standard_parallel: [-71.0]
straight_vertical_longitude_from_pole: [0.0]
version: 05-Nov-2019 (v1.38)
xmin: [-3333000]
ymax: [3333000]xarray.Dataset
- y: 13333
- x: 13333
- x(x)int32-3333000 -3332500 ... 3333000
- long_name :
- Cartesian x-coordinate
- standard_name :
- projection_x_coordinate
- units :
- meter
array([-3333000, -3332500, -3332000, ..., 3332000, 3332500, 3333000], dtype=int32) - y(y)int323333000 3332500 ... -3333000
- long_name :
- Cartesian y-coordinate
- standard_name :
- projection_y_coordinate
- units :
- meter
array([ 3333000, 3332500, 3332000, ..., -3332000, -3332500, -3333000], dtype=int32)
- bed(y, x)float32dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- grid_mapping :
- mapping
- long_name :
- bed topography
- source :
- IBCSO and Mathieu Morlighem
- standard_name :
- bedrock_altitude
- units :
- meters
Array Chunk Bytes 678.13 MiB 34.33 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type float32 numpy.ndarray - errbed(y, x)float32dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- grid_mapping :
- mapping
- long_name :
- bed topography/ice thickness error
- source :
- Mathieu Morlighem
- units :
- meters
Array Chunk Bytes 678.13 MiB 34.33 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type float32 numpy.ndarray - firn(y, x)float32dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- grid_mapping :
- mapping
- long_name :
- firn air content
- source :
- FAC model from Ligtenberg et al., 2011, forced by RACMO2.3p2 (van Wessem et al., 2018)
- units :
- meters
Array Chunk Bytes 678.13 MiB 34.33 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type float32 numpy.ndarray - geoid(y, x)int16dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- geoid :
- eigen-6c4 (Forste et al 2014)
- grid_mapping :
- mapping
- long_name :
- EIGEN-EC4 Geoid - WGS84 Ellipsoid difference
- standard_name :
- geoid_height_above_reference_ellipsoid
- units :
- meters
Array Chunk Bytes 339.07 MiB 17.17 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type int16 numpy.ndarray - mask(y, x)int8dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- flag_meanings :
- ocean ice_free_land grounded_ice floating_ice lake_vostok
- flag_values :
- [0, 1, 2, 3, 4]
- grid_mapping :
- mapping
- long_name :
- mask
- source :
- Antarctic Digital Database (rock outcrop) and Jeremie Mouginot (pers. comm., grounding lines)
- valid_range :
- [0, 4]
Array Chunk Bytes 169.53 MiB 8.58 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type int8 numpy.ndarray - source(y, x)int8dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- flag_meanings :
- REMA_IBCSO mass_conservation interpolation hydrostatic streamline_diffusion gravity seismic multibeam
- flag_values :
- [1, 2, 3, 4, 5, 6, 7, 10]
- grid_mapping :
- mapping
- long_name :
- data source
- source :
- Mathieu Morlighem
- valid_range :
- [1, 10]
Array Chunk Bytes 169.53 MiB 8.58 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type int8 numpy.ndarray - surface(y, x)float32dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- grid_mapping :
- mapping
- long_name :
- ice surface elevation
- source :
- REMA (Byrd Polar and Climate Research Center and the Polar Geospatial Center)
- standard_name :
- surface_altitude
- units :
- meters
Array Chunk Bytes 678.13 MiB 34.33 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type float32 numpy.ndarray - thickness(y, x)float32dask.array<chunksize=(3000, 3000), meta=np.ndarray>
- grid_mapping :
- mapping
- long_name :
- ice thickness
- source :
- Mathieu Morlighem
- standard_name :
- land_ice_thickness
- units :
- meters
Array Chunk Bytes 678.13 MiB 34.33 MiB Shape (13333, 13333) (3000, 3000) Count 26 Tasks 25 Chunks Type float32 numpy.ndarray
- Author :
- Mathieu Morlighem
- Conventions :
- CF-1.7
- Data_citation :
- Morlighem M. et al., (2019), Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nature Geoscience (accepted)
- Notes :
- Data processed at the Department of Earth System Science, University of California, Irvine
- Projection :
- Polar Stereographic South (71S,0E)
- Title :
- BedMachine Antarctica
- false_easting :
- [0.0]
- false_northing :
- [0.0]
- grid_mapping_name :
- polar_stereographic
- ice_density (kg m-3) :
- [917.0]
- inverse_flattening :
- [298.2794050428205]
- latitude_of_projection_origin :
- [-90.0]
- license :
- No restrictions on access or use
- no_data :
- [-9999.0]
- nx :
- [13333.0]
- ny :
- [13333.0]
- proj4 :
- +init=epsg:3031
- sea_water_density (kg m-3) :
- [1027.0]
- semi_major_axis :
- [6378273.0]
- spacing :
- [500]
- standard_parallel :
- [-71.0]
- straight_vertical_longitude_from_pole :
- [0.0]
- version :
- 05-Nov-2019 (v1.38)
- xmin :
- [-3333000]
- ymax :
- [3333000]
Compute the two components of the gradient using the numpy function
numpy.gradientThe second argument is the grid spacing, which is the same in both directions in this case. It is stored in the attributes of the xarray (and were loaded from the zarr in cloud storage).
dhdy, dhdx = np.gradient(bm.surface.values,bm.attrs['nx'][0])
4. Compute the slope by taking the pythagorean sum of the two components.¶
\[
S = \left(\left(\frac{\partial h}{\partial x}\right)^2 + \left(\frac{\partial h}{\partial xy}\right)\right)^\frac{1}{2}
\]
and add this to the xarray as a new field.
Also define a name and units for the new field for plotting and general documentation.
slope = xr.zeros_like(bm.thickness)
slope.values = np.sqrt(dhdy**2 + dhdx**2)
bm['slope'] = slope
bm.slope.attrs['long_name'] = 'ice surface slope'
bm.slope.attrs['units'] = '-'
5. Plotting¶
subset_of_slope = bm.slope.isel(x=slice(2000,3000), y = slice(4000,5000))
options = dict(x='x', y='y', clim=(0, 0.001), width=800, height=800, cmap='viridis')
subset_of_slope.hvplot(**options);
bm.slope.isel(x=slice(2000,6000), y = slice(4000,7000)).plot(vmin=0, vmax=0.005, size = 10);
all_slope_coarse = bm.slope.coarsen(x=5,y=5,boundary='trim').mean().persist()
options = dict(x='x', y='y', clim=(0, 0.001), width=800, height=800, cmap='viridis')
all_slope_coarse.hvplot(**options)